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ABSTRACT 

We study the formation of giant dense cloud complexes and of stars within them by means of 
SPH numerical simulations of the collision of gas streams ( "inflows" ) in the warm neutral medium 
(WNM) at moderately supersonic velocities. The collisions cause compression, cooling and turbulence 
generation in the gas, forming a cloud that then becomes self-gravitating and begins to collapse glob- 
ally. Simultaneously, the turbulent, nonlinear density fluctuations induce fast, local collapse events. 
The simulations show that: a) The clouds are not in a state of equilibrium. Instead, they undergo 
secular evolution. During the early stages of the evolution, their mass and gravitational energy E g 
increase steadily, while the turbulent energy reaches a plateau, b) When E g becomes compara- 
ble to £k, global collapse begins, causing a simultaneous increase in \E g \ and E^ that maintains a 
near-equipartition condition \E S \ ~ 2E^. c) Longer inflow durations delay the onset of global and 
local collapse, by maintaining a higher turbulent velocity dispersion in the cloud over longer times, 
d) The star formation rate is large from the beginning, without any period of slow and accelerating 
star formation, e) The column densities of the local star-forming clumps are very similar to reported 
values of the column density required for molecule formation, suggesting that locally molecular gas 
and star formation occur nearly simultaneously. The MC formation mechanism discussed here natu- 
rally explains the apparent "virialized" state of MCs and the ubiquitous presence of HI halos around 
them. Also, within their assumptions, our simulations support the scenario of rapid star formation 
after MCs are formed, although long (> 15 Myr) accumulation periods do occur during which the 
clouds build up their gravitational energy, and which are expected to be spent in the atomic phase. 

Subject headings: Instabilities — ISM: clouds — ISM: evolution — Shock waves — Stars: formation — 
Turbulence 



1. INTRODUCTION 

The evolution of molecular clouds (MCs) re- 
mains an unsolved problem to date. While tra- 
ditionally MCs have been thought of as virialized 
structures in the interstell ar med ium (I SM) (e.g., 
de Jong. Boland fc Dalgarno I Il980t IMcKee et al. I 
19931 Blitz fcWilliams 



2001 



IVazquez-Semadeni et al. I 



2006 



Ballesteros-Pare dcs & Hart mann I 2006; 
Ballesteros-Paredesl 120061 see also the reviews by 



garno lysu: McK.cc et al. I 
si 119991: IMcKee I 11999ft . with 



relatively long lifetimes ( Blitz fc Shu I [l980T ) and a 



significant delay before t hey b e gin forming stars (e.g. 
Palla fc Startler! l2000l 120021: iTassis fc Mouschovias I 



20041: iMouschovias. Tassis fc Kunz I 120061) . re- 
cent studies have suggested that MCs be- 
gin forming stars shortly after they them- 
selves form, and are non-equilibrium entities 
(Ballesteros-Paredes. Vazquez-Semadeni fc Scalo 1 119991 : 



Mac Low & Klessen 2004 and Ballesteros-Paredes et 
al. 2006). One approach that can shed light on this 
problem is to perform numerical simulations of the 
MC formation process by generic compressions in the 
warm neutral medium (WNM). Numerous studies of 
this kind exist, al though self-gravity has in general 
not been included dHennebelle fc Perault 1 119991 120001: 
20001 I2002t llnutsuka fc Koyama l 

120051 



Koyama & Inuts uka 
2004: I Audit fe Hennebeile 



2002; 
2005; 



Hcitsch et al. 



2006t IVazquez-Semadeni et al. 1 120061 hereafter Paper 



Ball esteros-Paredes. Hartmann fc Vazquez-Semadenil 



1999; 



Elmcgreen I 



Hartmann, Ballesteros-Paredcs & Bcrgin 
Vazquez-Semadeni, Ballesteros-Paredcs fc KlesserTI 
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I)- 

These studies have shown that a substantial fraction of 
the internal turbulence, with velocity dispersions of up 
to a few km s _1 , can be produced in shock-bounded lay- 
ers by a combined dynamical+ther mal instability, whose 
nature is still not well determined dHunter et al. 1 1 1986 ; 
I Walder fc Folini II1998L 120001: iKovama fc Inutsuka] 2002 : 
lAnmtBakan fc Kovamal 2004t [Heitsch et al. I 120051 . 12006 ; 
iFolini fc Walder 1 120061) . Paper I also showed that the 
dense gas, there defined as the gas with densities larger 
than 100 cm -3 , is at systematically higher thermal pres- 
sures than the mean interstellar value by factors of 1.5-5, 
implying that MC formation by colliding streams of dif- 
fuse gas can account at least partially for the clouds' 
excess pressure and turbulent nature. 

Only a few works have included self-gravity in sim- 
ulations of colliding flows. The pioneering work of 



2 



Vazquez-Semadeni et al. 



iHunter et al. I (|1986f ) considered the collision of su- 
personic streams within MCs, finding that the frag- 
ments produced by the dynamical instability (which they 
termed "Rayleigh- Taylor-like" ) were able to accrete mass 
until they became gravitationally unstable and collapse. 
However, this study was designed to investigate flows 
within MCs rather than the formation of the clouds 
themselves, and so it considered much smaller (sub- 
parsec) scales and different density, temperature and 
cooling regimes than are relevant for the process of MC 
formation. 

More recent numerical studies of MC formation 
not including self-g ravity ( Kovam a fe Inutsuka I l2002t 
iHeitsch et al. 1120051 ) have estimated that the individual 
clumps formed in the compressed layer formed by the 
colliding streams are themselves gravitationally stable, 
although the entire compressed-gas complex would be 
gravitationally unstable, had self-gravity been included. 

This raises the issue of whether the clumps can be 
considered in isolation. The question of which gas 
parcels are involved in the collapse of a particular 
object is a delicate one, and lies, for example, at the 
heart of the debate of whether competitive accretion 
is relevant or not in the process of star formation 
within molecular clouds: A model in which clumps are 
isolated dense objects immersed in a more tenuous, 
unbound medium leads to t he conclusion that competi- 
tive accretion is irrelevant (|Krumholz. McKee fc Klein I 
120051 ). Instead, com p etitiv e accretion plays a central 
role (IBonnell et al. I 119971: iKlessen. Bate fc Burkert I 
1998t iBonnell et all l2001t iSesscn & Burkert I ^"oq 



2001 



IKlessen I l2001allb F IBonnell fc Bate I 120061) in 
a scenario in which the local density peaks are 
the "tips of the iceberg" of the density distribu- 
tion, the flows are organized on much larger scales 
(Ballcsteros-Paredes, Vazqucz-Scmadcni & Scalo 1999; 
B alle steros-Paredes. Hartmann fc Vazquez-Semadenil 
1999tlBate. Bonnell fc Bromm II2003D . and the evolution 



and energetics of the density peaks cannot be considered 
in isolation, but as part of the global flow i nstead 
( Ballcsteros-Pared es. Vazquez-Semadeni fc Scalo I 
19991: 



Klcsscn, Hcitsch & Mac Low 



Heitsch. Mac Low fc Klessen I 



2000; 



2001 



Shadmehri. Vazquez-Semadeni fc Ballesteros-Paredesl 
20021 iTillev fc Pudritz 112004 120051 : iBallesteros -Paredes I 
20061 : iDib et all I2006D . In a similar fashion, then, 
accurate determination of whether gravitational collapse 
can be triggered by moderate-Mach number stream 
collisions in the WNM, requires to self-consistently 
include self-gravity in the numerical simulations. 

Simulations of molecular cloud formation by 
the passage of a self-gravitating , clumpy medium 
through Galactic spiral shocks ([Bonnell et al. I 120061 : 
iDobbs. Bonnell fc Pringle1l2006[ ) have shown that realis- 
tic velocity dispersions and densities can be generated in 
the resulting post-shock clouds, which are furthermore 
driven to produce local collapse events leading to star 
formation, although a large fraction of the mass in 
the clouds remains gravitationally unbound. However, 
because these simulations were isothermal, the gas 
had to be assumed to already be clumpy and at very 
low temperatures (~ 100 K) previous to the passage 
through the shock, while recent studies including 
cooling leading to thermal bistability of the atomic 



gas (iKovama fc Inutsuka I l2002t lAudit fc Hennebelle I 
120051 : IHeitsch et al. II2005L I2006L Paper I) suggest that 
the dumpiness is generated by the shock compression, 
which nonlinearly triggers the thermal instability, 
induces supersonic turbulence, and causes the formation 
of dense clumps in the compressed layers. 

It is thus essential to model the evolution of the cloud 
within the frame of its more diffuse environment includ- 
ing both self-gravity and cooling. Not only this is indis- 
pensible in order to model the generation of turbulence 
and dumpiness, but also to model the cloud's evolution 
as it exchanges mass and energy with its surroundings. 
Note that the generation of turbulence by the compres- 
sion itself implies that the turbulence within the cloud 
is driven while the large-scale converging motions per- 
sist, and decaying afterwards. This is a mixed regime 
that also requires self-consistent modeling of the cloud 
and its surroundings, rather than, for example, random 
Fourier driving in closed boxes. 

Finally, such a unified (or "holistic") description of 
mole cular cloud formation and evolution, as outlined 
by IVazquez-Semadeni. Ballesteros-Paredes fc Klessen I 
(|200a K including the star formation epoch, allows 
also the investigation of a number of key problems, 
such as the evolutionary and star-formation timescales, 
the evolution of the cloud's gravitational and kinetic 
energies, and whether it settles into virial equilibrium or 
not, and the effect of self-gravity on the cloud's physical 
conditions. 

In this paper we investigate numerically the evolution 
of interstellar gas as it collects from the WNM, shocks, 
suffers a phase transition to a cold and dense state, and 
finally begins to form stars, by means of numerical simu- 
lations of colliding gas streams in the WNM, in the pres- 
ence of self-gravity and cooling leading to thermal insta- 
bility. Since this is a problem that naturally collects large 
amounts of gas in small volumes, first by compression- 
triggered thermal instability and then by gravitational 
instability, it is convenient to use some Lagrangian nu- 
merical scheme, and we opt for a smoothed particle hy- 
drodynamics code, sacrificing the possibility of including 
magnetic fields. We also omit for now a self-consistent 
description of the stellar feedback on the gas and of the 
chemistry. Thus, our simulations will not adequately de- 
scribe the latest stages of MC evolution (in which stellar 
energy feedback may significantly influence the dynam- 
ics of the cloud) nor are they capable of predicting the 
transition from atomic to molecular gas accurately. We 
defer these tasks to subsequent papers. 

The plan of the paper is as follows. In £j2]we first give a 
qualitative description of the physical system, the main 
phenomena involved, and our expectations for the evolu- 
tion. In fj3]we describe the numerical code, the physical 
physical setting, and the choice of parameters. Then, in 
H] we describe the results concerning the evolution of the 
cloud's mass and its kinetic and gravitational energies, 
and the onset of star formation. In <j5] we discuss the 
implications and limitations of the recent work. Finally, 
in Sj5]we present a summary and draw some conclusions. 

2. QUALITATIVE DISCUSSION 

It is convenient to first present an overview of the sys- 
tem, its expected evolution, and the physical processes 
on which it relies (see also the description in Paper I). 
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Since we are interested in studying the formation of a 
dense cloud out of the diffuse medium, we take as initial 
condition a uniform-density region within the WNM, of a 
few hundred parsecs across, with an initially compressive 
velocity field consisting of two oppositely-directed inflows 
with speed v- m f and Mach number Ai ln {, with respect to 
the WNM temperature. This compressive field is gener- 
ically representative of either the general turbulence in 
the WNM, or of motions triggered by some large scale 
agent, such as gravitational or Parker instabilities, or 
shell collisio ns. Since the velocity dispersion in the WNM 
is transonic (jKulkarni fc Heiles Ifl987t iHeiles fc Trolaiidl 
2003j), we consider values of M. m f ~ 1. 

As is well known, the atomic ISM is thermally 
bistable, with two stable phases being able to coex- 
ist at the same pressure but different densities and 
temperatures, mediated by a thermally unstab le range 
(|Field 111965b iField. Goldsmith fc Habing lH96l) . The 
collision of transonic WNM streams produces a thick 
shock-bounded slab, in which the gas is out of ther- 
mal equilibrium between radiative heating and cool- 
ing. This shocked slab is nonlinearly thermally un- 
stabl e (jKovama fc Inutsuka 1 120001 : iKritsuk fc Norman"! 
2002), and as it flows towards the collision plane, 
it undergoes a phase transition to the cold neu- 
tral medium (C NM), causing the formation of a 
thin cold layer. dHennebelle fc Perault I I1999L 120001: 
Kovama fc Inutsuka 2000, 2003 llnutsuka fc Kovama I 
20041: lAudit fc Hennebeilel hooa IHeitsch et al. I 120051 
20061 Paper I). Subsequently, the boundary of this 
layer is destabilized, probably by a combination of 
thermal instabil ity and nonlinear thin-shell instability 
(|Vishniac II1994I ). producing fully-develop ed turbulence 
in th e layer for sufficiently large M\„f dHunter et al 

1 /~\ ci r\ 1 1 r _ _ _ o t .... j ... .1 _ I ITTTTTm It j _ i _ tr _ _ __ . _ 



19861 : JKovama fc Inutsuka 1 120021 : llnutsuka fc Kovama 
2004 IHeitsch et al. 1120051 12006L Paper I). As the layer 
becomes turbulent, it thickens, and becomes a fully 
three-dimensional structure (Paper I), to which we sim- 
ply refer as "the cloud". 

Moreover, the thermal pressure in the dense cloud is in 
pressure balance with the total (thermal + ram) pressure 
in the inflows, which is significantly larger than the mean 
thermal pressure in the ISM (Paper I). Consequently, the 
density in the cold layer can reach densities significantly 
larger than those typical of the CNM. Thus, this cloud 
is turbulent, dense, and overpressured with respect to 
the mean ISM pressure; i.e., it has properties typical of 
molecular clouds (Paper I) . The main missing link in this 
scenario is whether the atomic gas can be readily con- 
verted t o molecular, although several studies sugge st it is 
feasible (jBergin et al. 112004 : iGfover fc Mac Low 1 120061 ) . 
Since we do not incorporate any chemistry nor radiative 
transfer in our code, we cannot follow this transition self- 
consistently, and simply assume that once the dense gas 
has a sufficiently large column density (see below), it is 
rapidly converted to the molecular phase. 

In summary, the compression produces a turbulent 
cloud of dense and cold gas. As it becomes denser and 
colder, its Jeans mass decreases substantially. Moreover, 
its own mass is increasing, so it rapidly becomes much 
more massive than its own Jeans mass. Also, as the cloud 
becomes more massive, its total gravitational energy (in 
absolute value) |-E g | increases substantially, and eventu- 
ally overcomes even its turbulent kinetic energy, and the 



cloud begins collapsing as a whole. Because it contains 
many Jeans masses, local collapse events begin to occur 
at the nonlinear density peaks produced by the turbu- 
lence, if gravity locally overcomes all available forms of 
support. Because the density enhancements arc nonlin- 
ear, the local collapse events occur on shorter timescales 
than the global contraction. 

The magnetic field sho uld also be considered. 
iHennebelle fc Perault I (|2000f ) have shown that the tran- 
sition from diffuse to dense gas can occur even in 
the presence of magnetic fields aligned with or oblique 
to the direction of compression. Moreover, most evi- 
dence, both theoretical and observational, points towards 
MCs and giant molecular complexes being m agnetically 
super c ritical or marginally critical at least (|Crutcher~1 
19991: IHartmann. Ballesteros-Paredes fc Bergin 2001b 



Bourke et al. 112001b ICrutcher. Heiles fc Trolandl [2003) . 
Since the magnetically supercritical case is not qualita- 
tively different from the non-magnetic case, the above 
d iscussion is applicable to supercritical clou d s as w ell . 

IHartmann. Ballesteros-Paredes fc Bergin I (|2001l ) have 
approximately quantified the above scenario, estimating 
that the accumulation time from the WNM should be of 
the order of 10-20 Myr, the accumulated gas must come 
from distances up to a few hundred parsecs, and that at 
the end of the accumulation period, the gas should be 
becoming molecular and self-gravitating at roughly the 
same time, because the column density for self-shielding 
and formation of molecular gas, 



As. 



1 - 2 x 10 21 cm" 2 , 



(1) 



is very similar to that required for rendering the dense 
gas gravitationally unstable, 



A gi - 1.5 x 10 21 P 4 1/2 cnT 2 , 



(2) 

where P4 is the pres sure expressed in un its of 10 
K cm" 3 (see also iFranco fc Cox I 119861 ). That 
is, the gas is converted to the molecular phase 
roughly simultaneously with the onset of gravita- 
tional instability, explaining why non-star-forming 
MCs are rare, at least in the solar neighbo rhood 

(IHartmann . Balles teros-Paredes fc Bergin I l2001b 

IHartmann I l2003b iBallesteros-Paredes fc Hartmannl 
l2006h . Note, ho wever, that this near eq uality is a 
coincidence (e.g., lElmegreenl 11985b Il993al) . since A g i 
depends on the ambient pressure, and therefore non- 
star-forming molecular clouds may be more common 
in higher- pressure environments, possibly explaining 
observations that up to 30% of the G MCs appear 
to b e starless in some external galaxies (|Blitz et al. I 
2006). On the other hand, this result may be due in 
part to incompleteness effects. Indeed, recent Spitzer 
observations have uncovered young stellar objects in 
Galactic mol ecular cores previously believed t o be 
starless (e.g., ICraspi et al~ll2005b iRho et al. II2006D . 

In the remainder of the paper, we proceed to quantify 
this scenario by means of numerical simulations designed 
for this purpose. 

3. THE MODEL 

3.1. Parameters and numerical scheme 

We consider a cubic box of length Lbox on a side, ini- 
tially filled with WNM at a uniform density of no = 
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1 cm" 3 , with a mean molecular weight of 1.27, so that 
the mean mass density is po = 2.12 x 10 -24 g cm" 3 , which 
implies total masses of 6.58 x 10 4 M Q and 5.26 x 10 5 Af Q 
for the two box sizes used, Lb ox = 128 pc and 256 pc, 
respectively. The gas has a temperature To = 5000 
K, the thermal-equilibrium temperature at that density 
(cf. i)3.3|) . Within the box, we set up two cylindrical 
and oppositely-directed inflows, each of length £- ln {, ra- 
dius r inf = 32 pc, and speed v in { = ±9.20 km s~ , 
corresponding to a Mach number Ai[ n i = 1.22 with 
respect to the sound speed of the undisturbed WNM, 
cwnm = 7.536 km s" 1 (see fig. [IJ. For the two values of 
£i n f we consider here (48 pc and 112 pc), the masses of the 
cylindrical inflows are 1.13 x 10 4 M Q and 2.64 x 10 4 M , 
respectively. To trigger the dynamical instability of the 
compressed slab, a fluctuating velocity field (computed 
in Fourier space) is initially imposed throughout the box 
at wavenumbers 4 < fcLbox < 8, and having amplitude 
u rms .i. The streams are directed along the x direction 
and collide at half the x extension of the box. Table Q] 
gives the values of these parameters for the simulations 
presented here. The runs are named mnemonically as 
LxxxAvO. yy, where xxx indicates the physical box size in 
parsecs and O.yy indicates the value of w rms ,i in km s~ . 

We use the N-body+smo othed particle hydrodynam- 
ics (S PH) code GAD GET (iSpringel. Yoshida fc White"! 
1200 1| ). as modified by Uappsen et al. Y ( 20051 ) to include 
sink particles, a prescription to describe collapsed ob- 
jects without following their internal structure. We use 
1.64 x 10 6 particles for the runs with Lbox = 128 pc, and 
3.24 x 10 6 particles for the run with Tbox = 256 pc. Since 
the larger-box run contains 8 times more mass than the 
smaller-box ones but uses only twice as many particles, 
the minimum resolved mass in the larger-box run is four 
times larger than in the smaller-box run. Unfortunately, 
the computing resources available in our cluster did not 
allow us to further increase the number of particles. 

We take as units for the code a mass of IMq, a length 
£q = 1 pc and a velocity va — 7.362 km s _ , implying a 
time unit to = Iq/vq = 0.133 Myr. 

3.2. Sink particles 

Sink particles are created when a group of SPH parti- 
cles become involved in a local collapse event. At that 
point, the SPH particles are replaced by a sink particle 
whose internal structure is not resolved anymore. The 
sink particle inherits the total angular momentum and 
mass of the SPH particles it replaces, and henceforth 
moves only in response to its inertia and the gravitational 
force, although it can continue to accrete SPH particles, 
should they come sufficiently close to the sink particle. 

The creation of a sink particle requires a number of 
conditions to be satisfied. First, the local density should 
exceed some user-defined threshold value rithr- If this oc- 
curs, then the code computes the total mass and angular 
momentum of the group of particles above n t hr to deter- 
mine whether they are gravitationally bound. If they are, 
then they are collectively replaced by the sink. At sub- 
sequent times, further SPH particles can be accreted by 
the sink if they come within an "accretion radius" r acc of 
the sink and are gravitationally bound to it. In all cases 
we use nthr = 10 5 cm -3 and r acc = 0.04 pc. This value 
of r aC c amounts to roughly 1/3 of the Jeans length at 



7i t hr- Due to the limited mass resolution we can afford, 
the sinks should in general be considered star clusters 
rather than isolated stars, although some of them may 
actually correspond to single stars. Nevertheless, with 
this caveat in mind, we refer to the sinks as "stars" in a 
generic way. 

3.3. Heating and cooling 

The gas is subject to heating and cooling functions of 
the form 



T = 2.0 x 10~ 26 erg s 

A(T) , /-1.184xl0 5 
— — = 10 exp ' 



T + 1000 



+ 1.4 x 10" 2 V / Texp 



-92 
~~T 



cm 



(3) 



(4) 



These functions are fits to the va rious heating (r) and 
coolin g (A) processes considered by Kovama fc Inutsuka 



2000), as given by equation (4) of 
2002). 5 The resulting thermal-equilibrium 



Kovama fe Inutsuka 
pressure 

P eq , defined by the condition n 2 A = nT is shown in fig. [2] 
as a function of density. We have abandoned the piece- 
wise power-law fit we have used in all of our previous 
pape rs (e.g., Vazquez-S emadeni. Passot fc Pouauet 



1995 
1995 



Passot, Vazquez-Semadeni & Pouquet 



IVazquez-Semadeni. Passot fc Pououet I 19961: 
Vazquez-Semadeni. Gazol fc Scalol l2000t iGazol et al. I 
20011 : iGazol. Vazquez-Semadeni fc Kim I l2005f) . includ- 
ing Paper I, mainly because the present fit is applicable 
at densities typical of molecular gas. 

The usual procedure to apply cooling to the hydrody- 
namic evolution involves the consideration of the cooling 
rate in the Courant condition to restrict the simulation 
time step. For the problem at hand, the high densi- 
ties reached behind the shocks, for example, would im- 
ply an exceedingly small time step and the simulation 
would become computationally unfeasible. But this ap- 
pears completely unnecessary, because all this means is 
that the thermal evolution happens faster than the dy- 
namical one, and therefore, as far as the hydrodynamics 
is concerned, it is instantaneous. Therefore, it is more 
convenient to use an approximation to the thermal evo- 
lution of the gas, which constitutes an extension of that 
used in Paper I, and which allows us to simply correct 
the internal energy after the hydrodynamic step has been 
performed, with no need to adjust the timestep. 

Consider the thermal equilibrium temperature, T eq , as 
a function of the gas density, and the corresponding in- 
ternal energy density e eq (T eq is a well defined function 
of the density n as long as A(T) is monotonic. This is 
true within the range of temperatures occurring in the 
simulations). Consider also the time required to radiate 
the thermal energy excess (if T > T cq ; acquire the energy 
deficit, if T < T eq ) 



-'eq 



n 2 A — nT 



(5) 



5 Note that eq. (4) in Koyama & Inutsuka (2002) contains two 
typographical errors. The form used here incorporates the necesary 
corrections, kindly provided by H. Koyama. 
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Then, we compute the new internal energy density e', 
after a time step dt, as 



+ (e - e cq ) exp(-dt/r A ). 



(6) 



Note that, if the gas is cooling down (or heating up) 
rapidly, ta <§; dt, exp(— dt/rjs) — > and the gas immedi- 
ately reaches its equilibrium temperature, without ever 
overshooting beyond that value. Conversely, if the gas is 
at very low density or is close to the equilibrium temper- 
ature, ta 3> dt, and equation ([6]) becomes, 



e' = e- dt(n 2 k-nT), 



(7) 



where we, again, use the fact that A(T), as given by 
equation (j4|, is a monotonic function of the temperature. 

3.4. Gravitational instability in cooling media 

The standard Jeans instability analysis is modified in 
a medium in which the balance between heating and 
instantaneous cooling produces a net polytropic behav- 
ior, characterized by an effective polytropic exponent 7 e , 
which is the slope of the thermal-equilibrium pressure 
versus density curve (fig. [2|). In thi s case, the Jeans 
length is found to be given by (e.g.. lElmegreen IIT991I ; 
IVazquez-Semadeni. Passot fc PouquetT fT996) 



L, 



7 e 7T<r 



-1/2 



(8) 



where c is the adiabatic sound speed and 7 is the heat 
capacity ratio of the gas. Equivalently, this instability 
criterion can be described as if the effective sound speed 
in the system were given by c c g = \J j c kT/ /i, where k 
is Boltzmann's constant and /1 is the mean molecular 
weight. Note that, when j c — 0, the thermal pressure 
does not react to density changes and therefore is unable 
to oppose the collapse, rendering the medium gravita- 
tionally unstable at all scales (L j = 0) . In a more realis- 
tic situation where the cooling is not instantaneous, the 
Jeans length will be limited by the cooling length. In this 
work, the cooling has a finite characteristic timescale, 
given by eq. ([5|). 

At the initial uniform density no = 1 cm~ 3 and tem- 
perature T = 5000 K, 7e ~ and the Jeans length is also 
near zero. This can be seen in fig. [21 where the vertical 
dotted line shows the initial uniform density, at which 
the local tangent to the P eq vs. p curve is seen to have a 
nearly zero slope. However, when the compression pro- 
duces a dense cloud at the midplane of the box, the rest 
of the gas decreases its density and enters the fully stable 
warm phase. For reference, at a density n — 0.5 cm~ 3 



and a corresponding equilibrium temperature T e , 



6300 



K, 7o = 0.834 and Lj = 1324 pc, implying a Jeans mass 
Mj = 3.47 x 10 7 M Q . Under these conditions, our runs 
with Lbox = 128 and Lbox = 256 pc would initially con- 
tain ~ 10~ 3 and ~ 10~ 2 Jeans masses, respectively. 

4. RESULTS 

4.1. Mass evolution 

In general, the simulations evolve as described in Pa- 
per I and in sec. $2] As an illustration, figs. [3] and |4] 
show selected snapshots of run L256Av0.17 viewed edge- 
on and face-on, respectively. In the electronic version of 



the paper, these figures are also available as full-length 
animations. The times indicated in the frames are in the 
code's internal unit, to — 0.133 Myr. 

The animation shows that the gas at the collision site 
begins to undergo a phase transition to the cold phase 
(the CNM) at roughly one coo l ing time after the collision 
(iHennebelle &: Perault 1 119991 ; IVazquez-Semadeni et al. I 
|2006[ ). becoming much denser and colder. Because this 
gas is in equilibrium with the total (thermal+ram) pres- 
sure of the inflow, the density overshoots far beyond that 
of standard CNM, well into the realm of molecular gas, 
with mean densities of several hundred cm -3 and tem- 
peratures of a few tens of Kelvins. This implies a much 
lower Jeans mass than that in the conditions of the ini- 
tial WNM and quickly the cloud's mass exceeds its Jeans 
mass. 

Defining the "cloud" as the material with n > 
50 cm -3 , the cloud's mass becomes larger than its mean 
Jeans mass at t ~ 2.2 Myr in runs with Lbox = 128 
pc and at t ~ 3 Myr in the run with Lbox = 256 pc, 
at which times the mean density in the dense gas is 
(n) ~ 65 cm -3 , the mean temperature is (T) ~ 50 K, 
the Jeans mass is Mj ~ 1200 M , and the free-fall time 
is Tff = Lj/c ~ 11 Myr, where c is the local sound speed. 
For comparison, the total mass in the two cylinders, 
which contain the material that initially builds up the 
cloud, is 1.13 x 10 4 M Q in the runs with Lbox = 128 pc, 
and 2.26 x 10 4 M Q in the run with Lb ox = 256 pc. Thus, 
the cloud eventually becomes much more massive than 
its mean Jeans mass. For example, in run L256Av0.17 
at t = 17.26 Myr, the (mass-weighted) mean density of 
the dense gas is (n) ~ 600 cm' 3 , the mean temperature 
is (T) ~ 37 K, and 7 C ~ 0.75. These values imply a 
mean effective Jeans mass Mj ~ 300 M Q , while at that 
time the cloud contains ~ 2.35 x 10 4 M@. Note that 
the cloud+sink mass might exceed the mass in the in- 
flows, since ambient material surrounding the cylinder 
is dragged along with the flow, effectively increasing the 
amount of gas in the inflow (see below). 

The animations show that global collapse does indeed 
occur, although at significantly later times than when 
the cloud's mass becomes larger than the Jeans mass, 
probably reflecting the role of turbulence as an opposing 
agent to the collapse. Indeed, the clouds in the 128-pc 
boxes begin global contraction at t ~ 8 Myr, while the 
cloud in the 256-pc box, with more than twice the inflow 
duration, begins contracting at t ~ 12.5 Myr. 

Figure[5] (top panel) shows the evolution of the cloud's 
mass in run L256Ai;0.17, together with the mass in col- 
lapsed objects (sinks), to be discussed further in §4.51 
It is clearly seen that the cloud's mass is not constant, 
but rather evolves in time. Gas at cloud densities first 
appears after 2 Myr of evolution, and the cloud's mass 
increases monotonically until t ~ 19.5 Myr, at which 
time it reaches a maximum value M cl - 2.7 x 1O 4 M . 
After this time, it begins to decrease because of the rapid 
conversion of gas to stars. 

By the end of the simulation, the total mass in stars 
plus dense gas is almost three times the initial mass in 
the inflows. In practice, the inflows last longer than the 
simple estimate Tj n f = Anf/finf would indicate, and in- 
volve more mass than the mass initially within the cylin- 
ders. This is due to the fact that, as the cylinders be- 
gin advancing, they leave large voids behind them that 
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have the double effect of slowing down the tails of the 
inflows and of dragging the surrounding gas behind the 
cylinders. For v- m i — 9.20 km s _1 , Ti n £ = 5.22 Myr for 
runs with £ inS = 48 pc (L128Aw0.24 and L128Av0.66), 
and Ti n f = 12.17 Myr for the run with ^; n f = 112 pc 
(L256Av0.17). Instead, in the first two runs, the inflow 
actually lasts ~ 11 Myr, while in the latter run it lasts 
~ 22.5 Myr. This can be observed in the animation of 
fig. [31 in which the velocity field arrows clearly show the 
gradual decay of the inflows, and their longer duration 
compared to the simple linear estimate. 

4.2. Cloud disruption 

Our simulations do not include any feedback 
from the stellar objects, while this process is prob- 
ably essential for the energy balance and pos- 
sibly the dest ruction of the cloud (e.g., ICox I 



19831 119851: iFranco. Shore fc Tenorio-Tagle I 



1994 



Hartmann, Ballestcros-Paredes & Bcrgin 2001 
Matzner 1 120021: INakamura fc Li I 120051 : iLi fc Nakamura I 
20061: iMellema et al. I l2006[ ). Thus, the simulations 
cannot be safely considered realistic after their stellar 
content would likely disrupt the cloud. We now present 
an a posteriori estimation of the time at which this 

o ccurs. 

IFranco. Shore fc Tenorio-Tagle I (|1994f ) have suggested 
that the maximum number of OB stars that the cloud can 
support at any time is given by the number of H II regions 
required to completely ionize the cloud. In particular, 
they found that the maximum number of massive stars 
that can be formed within a molecular cloud is 



N, 



OB, inside 



l6M cA nf 7 



(9) 



where M Cj 4 is the mass of the cloud, in units of 10 4 
Mq, 77,3 is the number density in units of 10 3 cm~ 3 , 
i<48 is the photoionizing flux from the massive stars, in 
units of 10 48 s -1 , c S! i5 is the isothermal sound speed 
in units of 15 km s _1 , and £ms,7 is a characteristic OB 
star ma in sequence lifetime, in units of 10 7 y r. How- 
ever, as IFranco. Shore fc Tenorio-Tagle I (|1994h pointed 
out, clouds arc more efficiently destroyed by stars at 
the cloud edge because the lower external pressure en- 
sures that the ionized gas expands rapidly away from 
the cloud, driving a fast ionization front into the dense 
material. Thus, in this case, these authors find that the 
maximum number of OB stars that the molecular cloud 
can support is given by: 



Nqb 



edge 



3M C , 4 n\ r ° 



FU 5 (c s>15 <MS,7) 6 / 5 



(10) 



In order to determine the time at which the MC in 
our simulation has formed a sufficient number of stars to 
disrupt it, we proceed as follows. Since the individual 
sink particles in the simulations in general correspond to 
clusters rather than to single stars, we cannot use the 
sink particle masses directly. Instead, for each temporal 
output of our simualtion , we fit a standard initial mass 
function (|Kroupa 1120011 ) to the total mass in sink par- 
ticles at that time, to obtain the the "real" distribution 



of stellar masses i n the simulation at that time. 6 W e 
then use Table 1 of lDi'az-Miller. Franco fc Shore I (| 19981) . 
which gives the ionizing flux as a function of stellar mass, 
to estimate the flux associated to each mass bin, and then 
integrate over all relevant masses to obtain the total F48 . 

Now, from eqs. or (fit)]) one can solve for M C) 4 as a 
function of Nob, to find the minimum cloud mass that 
survives the ionizing radiation from the existing OB stars 
in the simulation. When this minimum surviving mass 
is larger than the actual cloud mass, we expect that the 
cloud will be disrupted. Since we define the "cloud" as 
the gas with densities above 50 cm -3 , we take n% = 0.05. 
Also, for simplicity we assume c S) is £ms,7 ~ 1- 

For run L256Aw0.17, fig. [H] shows the evolution of the 
cloud's mass together with the minimum disrupted cloud 
mass under the two estimates, eqs. ([9]) and (fTQ|) . We see 
that the cloud is expected to be disrupted at t = 20.3 
with either estimate. 
At this time the mass in sinks is M s i n k s = 5.25 x 10 3 M m 



Zsinks = iu"m 

d = 2.67 x 10 4 Mq, - " ■ ' 
star formation efficiency (SFE) up to this time has been 



and the cloud's mass is M r \ = 2.67 x 10 4 Mp i , so that the 



SFE = 



16%. 



(11) 



This number is still l arger than observat ional estimates 
for cloud complexes (|Mvers et al. I IT986) . and suggests 
that additional physical processes, such as longer inflow 
durations ( §4.6p or magnetic fields ( §5.2|) may be neces- 
sary to further reduce the SFE. 

We conclude from this section that after t rs 20.3 Myr 
run L256Ai>0.17 is probably not realistic anymore, as far 
as the evolution of a real MC is concerned, and we thus 
restrict most of our subsequent discussions to times ear- 
lier than that, except when later times may be illustrative 
of some physical process. Self-consistent inclusion of stel- 
lar fe^db^ck_Jn_the simula tions , similarly to the studie s 
of INakamura fc Li I (|2005h and ILi fc Nakamura I (|2006l) . 
to investigate the final stages of MC evolution will be the 
subject of future work. 

4.3. Turbulence evolution 

The collision also generates turbulence in the dense 
gas, and so the cloud can be considered to be in a driven 
turbulent regime while the inflows persist, and in a decay- 
ing regime after the inflows subside, although the tran- 
sition from driven to decaying is smooth, rather than 
abrupt, since the inflows subside gradually, as explained 
above. Figure [5] (bottom panel) shows the evolution of 
the velocity dispersion a for the dense gas along each of 
the three coordinate axes for run L256At>0.17. During 
the first 10 Myr of the cloud's existence (2 < t < 12 Myr) 
o~y ~ c z ~ 4 km s _ , although it is interesting to note 
that the turbulence appears to be anisotropic, with o~ x 
first increasing and then decreasing. It reaches a maxi- 
mum of more than twice as large (~ 8.5 km s _1 ) as that 
of u y and a z at t w 6 Myr, reflecting the fact that the 
inflows are directed along this direction. This suggests 
that the generation of transverse turbulent motions is not 



6 Note that this estimate for the number of massive stars as a 
function of time is probably a lower limit, as recent simulations 
suggest that massive stars in clusters tend to form first, while stars 
formed later tend to be less massive because o f the competitio n for 
accretion with the existing massive stars (e.g. [Bonncll 200l). 
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100% efficient. Of course, this effect is probably exag- 
gerated by our choice of perfectly anti-parallel colliding 
streams. Clouds formed by obliquely colliding streams 
are likely to have more isotropic levels of turbulence. 

After t ~ 6 Myr, a x begins to decrease, reflecting 
the weakening of the inflows, settling at o~ x ~ 6 km s _1 
by t ~ 12 Myr. At this time, o~ v and o~ z begin to in- 
crease, reflecting the global collapse of the cloud on the 
yz plane, while o~ x remains nearly stationary, unaffected 
by the global planar collapse, until t ~ 18 Myr, at which 
time it also begins to increase. This coincides with the 
onset of star formation, and probably reflects the local, 
small-scale, isotropic collapse events forming individual 
collapsed objects. 

4.4. Evolution of the density, pressure and temperature 
distributions 

An important consideration for understanding the pro- 
duction of local collapse events is the distribution of the 
density and temperature as the cloud evolves. This is 
shown in fig. [7) where the normalized mass- weighted his- 
tograms of density (left panel) and of temperature (right 
panel) are shown at times 4, 8, 12, 14, 16, and 20 Myr 
for run L256Av0.17. It can be seen that a cold phase 
already exists at t = 4 Myr, although at this time the 
cloud's conditions do not greatly exceed those of the stan- 
dard CNM, meaning that turbulence is only building up 
in the cloud at this stage, at which the simulation re- 
sembles a two-phase medium. At t = 8 Myr, the density 
maxima and temperature minima have shifted to more 
extreme values, which persist up to t — 12 Myr, indi- 
cating that the density and temperature fluctuations are 
predomin antly created by the supersonic turbulence in 
the cloud (|Vazquez-Semadeni et al. I f2006). However, by 
t = 14 Myr, the density and temperature extrema are 
seen to be again moving towards more extreme values, 
and this trend persists throughout the rest of the sim- 
ulation, indicating that gravitational collapse has taken 
over the evolution of the density fluctuations (see also 
SU). 

It thus appears that in this particular simulation, the 
turbulent density fluctuations act as simply as seeds 
for the subsequent growth of the fluc tuations by grav- 
itation al instability, as proposed by IClark fc Bonnelil 
(2005). However, both in that paper as in the present 
study, the global turbulence is already decaying, and by 
definition it cannot then continue competing with self- 
gravity. Minimally, a fundamental role of turbulence, 
even in this decaying state, must be to create nonlin- 
ear density fluctuations, of much larger amplitudes than 
would be created by thermal instability alone, which can 
collapse in shorter times than the whole cloud (cf. ^4.51 
and 15. ip . In any case, simulations in which the inflow 
lasts beyond the onset of collapse would be desirable, but 
unfortunately, as mentioned in i)4.61 we have not found 
it feasible to attempt them with the present code and 
computational resources, and such a study must await a 
different numerical scheme and physical setup. 

4.5. Energy evolution and star formation 

Figure [5] (top panel) shows the evolution of the ki- 
netic (i?k) and thermal (-Eth) energies within a cylinder 
of length 16 pc and radius 32 pc centered at the middle 



point of the numerical box, together with the evolution of 
the (absolute value of the) gravitational energy (\E g \) for 
the entire simulation box. This cylinder contains most 
of the cloud's mass thoughout the simulation, although 
it also contains sizable amounts of interspersed WNM. 
Figured] (bottom panel), on the other hand, shows the 
evolution of the same energies but with E^ and -Etii cal- 
culated for the dense gas only. 

The gravitational energy is shown for the entire box 
because of the practical difficulty to evaluate it only for 
the cloud's mass, although we expect the latter to domi- 
nate the global gravitational energy when the cloud has 
become very massive. Indeed, we observe that |J5 g | at 
late times is very large and dominates the other forms 
of energy. A brief period of positive values of \E g \ is ob- 
served for < t < 9 Myr, which can be understood as 
follows: Because the boundary conditions are periodic, 
Poisson's equation for the gravitational potential is ac- 
tually modified to have the density fluctuation p — (p) as 
its source. This is standard fare in cosmological simula- 
tions (see, e.g., lAlecian fe Leoratl[l9 88; Kolb fc Turner I 
1990). Moreover, gravity is a long-range force, and so 
large-scale features tend to dominate E g . Thus, during 
the early epochs of the simulation, when the cloud is be- 
ginning to be assembled, the dominant density features 
are the voids left behind by the cylinders because they 
are large, while the cloud is small and not very massive. 
However, once the cloud becomes sufficiently massive, 
and the voids have been smoothed out, the cloud domi- 
nates E g . The period of positive E g is omitted from both 
panels of fig. [5] 

Comparison of the energy plots (figs. [8] top and bot- 
tom) with the evolution of the gas and sink masses (fig. 
[5] top) and the animations shows some very important 
points. First and foremost, the cloud is never in virial 
equilibrium over its entire evolution. Instead, we can 
identify three main stages of evolution. First, over the 
period 8.5 < t < 14 Myr, \E g \ increases monotonically, 
transiting from being negligible compared to E^ and E t h 
to becoming larger than either one of them. The exact 
time at which this occurs depends on what system is be- 
ing considered. It occurs at t ~ 10 Myr when only the 
energies in the dense gas are considered, while it occurs 
at t ~ 18 Myr when the entire cylinder, which includes 
substantial amounts of warm gas, is considered. Since 
global gravitational contraction starts at t ~ 12 Myr, it 
appears that the true balance is bracketed by the esti- 
mates based on the full cylinder volume and on the dense 
gas, though closer to the latter. Thus, over the interval 
8-5 < t < 12 Myr, the increase in \E g \ is driven mostly 
by the cooling and compression of the gas. 

Second, from t ~ 12 Myr to t ~ 24 Myr, \E g \ continues 
to increase (E g becomes more negative), but now driven 
by the global collapse of the cloud. Over this period, we 
see that E^ for the dense gas (fig. bottom panel) closely 
follows \E g \, approximately satisfying the condition 

\E g \ = 2E k . (12) 

We see that in this case, this condition is a signature of 
collapse, not of gravitational equilibrium, even though it 
lasts for nearly 12 Myr. The maintenance of this equipar- 
tition arises from the fact that the cloud is converting 
gravitational potential energy into kinetic energy, but is 
replenishing the former by the collapse itself. 
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Note that we have cautioned in M.2\ that our simula- 
tions may not be realistic after t ~ 20.3 Myr, because 
by that time the stellar energy feedback may be suffi- 
cient to revert the global collapse and disperse the cloud, 
but even in that case the cloud will have evolved out of 
equilibrium up to that point. 

Finally, after t ~ 24 Myr, the gravitational and kinetic 
turbulent energies saturate and begin to vary in approx- 
imate synchronicity. Although this late stage may not 
be representative of actual MCs, it is important to un- 
derstand what is happening in the simulation. The near 
constancy of \E g \ and E^ could naively be interpreted 
as a final state of near virial equilibrium. However, in- 
spection of the animations shows otherwise. At t ~ 24 
Myr, the main body of the cloud is completing its global 
collapse, and the mass in stars is beginning to exceed the 
mass in dense gas, which is itself decreasing. So, at this 
point, the gravitational energy is beginning to be domi- 
nated by the stars, rather than by the gas. Moreover, the 
face-on animation (fig. 2} shows that at this time (t = 180 
in code units, shown in the animations) the star cluster is 
beginning to re-expand, after having reached maximum 
compression. This leads to the decrease of \E S \ between 
t ~ 24 and t ~ 28 Myr. Meanwhile, E^, which is domi- 
nated by the dense gas, decreases because the dense gas 
is being exhausted. Nevertheless, the outer "chaff" of 
the dense gas, mostly in the form of radial filaments, is 
continuing to fall onto the collapsed cloud. This leads to 
a new increase in both energies, as this is a secondary 
collapse. This situation repeats itself at t ~ 32 Myr 
(t = 244 in code units) after the secondary collapse ends 
and its second-generation star cluster begins to expand 
away. The end of each collapse and the re-expansion of 
the clusters is marked by kinks in the graph of Af s ; n k s 
versus time in fig. [5] {top), indicating a decrease in the 
star formation rate. 

It is important to recall as well that the condition 
given by eq. (|12[) is not sufficient for virial equilibrium. 
The necessary and sufficient condition for this is that 
the second time derivative of the moment of inertia of 
the cloud be zero, and eq. (I12p cannot guarantee this, 
as ma ny other terms enter the full virial balance of the 
cloud (fBallesteros-Paredes. Vazauez -Scmadc ni fc Scalo I 
ll999tlBallesteros-Paredesll2006tlDib et al 112006ft . Thus, 
the observed closeness of \E g \ and E^ in actual molecu- 
lar clouds must exclusively be considered an indication 
of near equipartition and pr obably of collapse, b ut not 
of virial equilibrium (see also iKlessen et al. 1 12005) . 

Thus, rather than evolving in near virial equilibrium, 
the cloud evolves far from equilibrium. \E S \ starts out 
from essentially zero and increases monotonically until it 
catches up with the thermal and kinetic energies. From 
that time on, gravity dominates the energy balance, lead- 
ing to collapse, which in turn causes a near equipartition 
between \E g \ and J^k, although both energies continue to 
vary systematically. 

Finally, figure [5] (top) also shows that star formation 
begins at t ~ 17.2 Myr, roughly 5 Myr after global con- 
traction of the cloud has begun. Yet, the star-forming 
local collapse events proceed much more rapidly than 
the global collapse of the cloud. This indicates that the 
turbulence in the cloud has created nonlinear density 
fluctuations, whose local free-fall time is shorter than 
that of the whole cloud. In particular, the mass in stars 



increases from zero to ~ 15% of the cloud's mass in ~ 3 
Myr (from t = 17 to t = 20 Myr). 

The main conclusions from this section are that a) the 
cloud evolves far from equilibrium all the way from its 
inception through its final collapse; b) as soon as E g 
dominates the dynamics, the cloud develops equiparti- 
tion indicative of the collapse, not of equilibrium, and c) 
star formation is rapid in comparison with the evolution 
and collapse of the whole cloud. 

4.6. Effect of the inflow duration and initial velocity 
dispersion 

In previous papers (IVazquez-Semadeni. Kim fe Ballesteros-Parede 



[2005T IVIzTni" ez-Semadeni et al. 1 12006ft . it has been ar- 
gued that the star formation efficiency is a sensitive 
function of whether the turbulence in the cloud is in a 
driven or in a decaying regime. Our system is driven 
at first, and gradually transits to a decaying regime, as 
pointed out in M4.3I In all of the runs in this paper, 
global cloud contraction and the subsequent star forma- 
tion occur after the inflows have weakened substantially 
(t > Ti n {). The main motivation behind run L256Av0.17 
was precisely to model an inflow of as long a duration 
as possible, to approximate the case of a driven cloud, 
althou gh this goal was not completely achieved. For ex - 
ample. lHartmann. Ballesteros-Paredes fc Bergin I (|2001l ) 
suggest accumulation lengths of up to 400 pc. Runs 
with even larger boxes (e.g., 512 pc) would be desirable, 
but they are either prohibitely expensive, or have an 
excessively poor mass resolution. In addition, vertical 
stratification effects would have to be considered. This 
will require a transition to the code Gadget2, which 
allows for non-cubic boundary conditions and/or inflow 
boundary conditions, a task we defer to a subsequent 
paper. 

Nevertheless, comparison of the runs with different in- 
flow lengths does shed light on the effect of a longer driv- 
ing duration. The cloud in run L256Aw0.17, whose inflow 
has £- m f = 112 pc and lasts ~ 22 Myr (although it begins 
weakening at T ~ r m t — 12.2 Myr), begins contracting 
at t ~ 12.5 Myr and begins forming stars at t ~ 17 Myr, 
while both runs with f; n f = 48 pc, for which the inflow 
lasts 11 Myr (begins weakening at t ~ r in f = 5.2 Myr), 
start to contract at t ~ 8 Myr (recall that the inflows 
have the same density and velocity in all runs). Run 
L128Av0.24 begins forming stars at t = 15.75 Myr, and 
run L128A?;0.66 does so at t = 13.77 Myr. 

These results clearly suggest that a longer inflow du- 
ration delays the onset of both global collapse and star 
formation, in spite of the fact that the cloud formed by 
it is more massive. This is attributable to the fact that 
the turbulence in the cloud is continuously driven by the 
inflow. This is verified by comparing the three compo- 
nents of the velocity dispersion for runs L256Aw0.17 and 
L128Aw0.24, respectively shown in fig. \5\ (bottom) and 
fig. El Both runs show an initial transient peak in a x 
ending at t rs 4 Myr. However, after this transient, <j x 
in run L256Au0.17 increases again, reaching a maximum 
of a x ~ 8.5 km s _ , and decreasing afterwards, until it 
nearly stabilizes at a value a x ~ 6 km s~ x at t ~ 11.5 
Myr. During this time interval, a y ~ a z ~ 4 km s~ . 
Instead, in run L128Aw0.24, a x does not increase again 
after the initial transient, and remains at a much more 
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moderate level of o~ x ~ 5 km s , while o~ y ~ o~ z ~ 
3 km s _1 . Clearly, the turbulence level is significantly 
higher in the longer-inflow run. 

Note that the delay in the onset of global collapse and 
of local star formation cannot be attributed to the ampli- 
tude of the initial velocity fluctuations in the inflow, since 
runs L128Av0.24 and L256Av0.17, which have compara- 
ble amplitudes, differ substantially in these times, while 
runs L128Au0.24 and L128Au0.66 begin global contrac- 
tion at almost the same time and differ by only ~ 15% 
in the time at which they begin forming stars, in spite 
of one having more than twice the velocity fluctuation 
amplitude of the other. This leaves the inflow duration 
as the sole cause of delay of both the global and local 
collapses. 

4.7. Local column density and star formation 

The initial conditions in our simulations are assumed 
to consist exclusively of atomic gas. However, the final 
density and temperature conditions are typical of molec- 
ular clouds, so molecule formation must occur somewhere 
along the way in the evolution of the clouds. One impor- 
tant shortcoming of our simulations is that they can- 
not distinguish between atomic and molecular gas, as no 
chemistry is included. Thus, we cannot directly tell from 
the simulations how long does it take for star formation 
to begin after molecular gas forms. We can, neverthe- 
less, measure the column density of the dense gas in 
star-forming regions, and c ompare it to t h e esti mates 
for self-shielding given by iFranco fe Cox I (119861) and 
lHartmann. Ballesteros-Paredes fc Bergin I (|2001l ). under 
solar Galactocentric conditions of the background UV 
field. The former authors find a threshold column 
desity for self-shielding of N ss ~ 5 x 10 20 cm~ 2 , 
while the latter authors quote values N ss ~ 1- 
2xl0 21 cm" 2 from Ivan Dishoeck fe Black I (|1988ft and 



Ivan Dishoeck fe Blake I ( 19981 ). So, we take a reference 



value of N st 



cm 



For comparison, in Table [2] we report the column den- 
sities of the dense gas in the first four regions of star 
formation in run L256Au0.17 at the time immediately 
before they begin forming stars. We see that the column 
densities fall in the range 0.5-2 xlO 21 cm -2 , suggesting 
that star formation locally occurs nearly simultaneously 
with the conversion of the gas from atomic to molecular. 

Moreover, we can compare the time at which we expect 
the bulk of the cloud to become molecular (i.e., the time 
at which the mean column density in the cloud is reaching 
N ss ) with the time at which star formation begins in the 
cloud. A lower limit for the time at which the mean 
column density in the cloud equals -/V ss is given by 

tss > = ^ M yr, (13) 

~ 2n v ini 5.8 x 10 19 cm" 2 y ' V ' 

where n$ and i?i n f have been defined in i j3.ll This es- 
timate is a lower limit because, as mentioned in i)4.11 
t>i n f is not constant, but actually decreases in time. For 
N ss = 10 21 cm -2 , we obtain t ss > 17.3 Myr, suggesting 
that the bulk of the cloud is still expected to be predom- 
inantly atomic by the time star formation is beginning. 

Together, these simple estimates suggest that stars 
form roughly simultaneously with the conversion of gas 
from atomic to molecular, with the local star form- 
ing regions being more advanced in the conversion pro- 



cess than the bulk of the cloud. In turn, this suggests 
that our model cloud would be observed as a collec- 
tion of predominantly molecular clumps immersed in a 
large, predominantly atomic, substrate. This sugges- 
tion is consistent with recent observations that substan- 
tial amounts of atomic gas coexist w i th the molecular 
phase in MCs dLi fe Goldsmith I [2001 iGoldsmith fe Li I 
l2005t iKlaassen et al. I I2005D. with perhaps even warm 
gas mixed in (|Hennebelle fc Inutsuka Il2006l ). Of course, 
a more conclusive confirmation of these estimates must 
await simulations in which the chemistry and radiative 
transfer are properly taken care of. 

5. DISCUSSION 

5.1. Implications 

5.1.1. Dynamic evolution and global versus local collapse 

The results presented in <2] pr ovide genera l suppor t 
to the scenario outlined in lElmegreen I (Tl993rj|) ; 
iBallesteros-Paredes. Vazquez-Semadeni fc Scalo I 



1999); Bal lesteros-Paredes. Hartmann fc Vazquez-Semadenil 



Hartmann, Ballesteros-Paredes & Bergin 



2001); Vazquez-Semadeni, Ballesteros-Paredes & Klcssen 
(120031); lHartmann I (120 03): Ma c Low fc Klessen I (120041); 
iHeitsch et al. I (120051. I2006D; I Vazquez-Semadeni et al. I 
(120061 ): IBallesteros-Paredes et al. I (|2006l) . The clouds in 
our simulations are never in a state of virial equilibrium 
before they convert most of their mass into stars. In- 
stead, they are in a continuously evolving state, initially 
obtaining their mass and turbulence simultaneously as 
they form out of the compression and cooling of diffuse 
warm gas (see also iHennebelle fc PerauTtl I1999L [2000; 
iKovama fc Inutsuka I l2002t llnutsuka fc KovamaTl2004l) . 
The gas initially has negligible self-gravity compared to 
its thermal support, but it quickly becomes super- Jeans 
because of the compression, the cooling, and the mass 
increase of the dense gas, until its self-gravitating 
energy eventually becomes comparable with the sum 
of its thermal and turbulent energies, at which point 
it begins to undergo global collapse. After this time, 
gravity becomes the main driver of the large-scale 
motions in these (semi-) decaying simulations, causing 
a near-equipartition, |_B g | ~ 2£"k, which however is 
indicative of collapse, not equilibrium. 

The nonlinear density fluctuations induced by the tur- 
bulence collapse faster than the whole cloud, as they have 
shorter free-fall times, and star formation proceeds vig- 
orously before global collapse is completed. This result 
is in contrast with the standard notion that linear den- 
sity fluctuations cannot lead to fragmentation because 
the fastest growing mode of gravitational instability in a 
nearly uniform medium is an overall contraction of the 
whole medium (|Larson I Il985f) . Thus, a crucial role of 
the turbulence in the medium is to create nonlinear den- 
sity fluctuations that have shorter free-fall times than the 
entire cloud. 

These results bring ba ck the scenario of globa l 
cloud collapse proposed bv iGoldreich fc Kwan I (|1974h . 
but with a twist t hat avoids the criticism of 
IZuckerman fc Palmer] (|1974[ ). namely that MCs could 
not be in global collapse because the star formation rate 
would be exceedingly high. Our simulations suggest that 
MCs may be undergoing global gravitational contraction, 
but the efficiency is reduced because local collapse events, 
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which involve only a fraction of the total cloud mass, pro- 
ceed faster than the global collapse. Once a sufficiently 
large number of stars have formed, their energy feedback 
may partially or completely halt the collapse. Further re- 
duction may be provided by supercritical magnet ic fields 
(IVazouez-Semadeni. Kim fc Ballesteros-Paredes I 120051 : 
iNakamura fc Li Il2005ft~ This scenario is consistent with 
the recent proposal of lHartmann fc Burkert I (|2006[ ) that 
the Orion A cloud may be undergoing gravitational col- 
lapse on large scales. 

5.1.2. Cloud's mass variation and cloud boundaries 

The fact that the cloud's mass is not constant is equiv- 
alent to the property that the locus of a Lagrangian 
boundary of the cloud defined (by means of a threshold 
density) at any given time does not remain at the cloud's 
boundary as time progresses, but instead it is later in- 
corporated into the interior of the cloud. This implies 
that the cloud does not always consist of the same set 
of fluid parcels, and that in a virial balance analysis of 
the cloud, there is non-zero flux of the physical variables 
through Eulerian boundaries defined at any given time 
(IBallesteros-Paredes. Vazouez-Semadeni fc Scalo 1 119991 : 
Shadmchri, Vazq uez-Semadeni fc Ballcst eros-Paredea 
2002t iDib et al Il2006t IBallesteros-Paredes 112006ft ~ 



5.1.3. Absence of slow, accelerating star formation phase 

It is important to note that the star formation rate 
(given by the slope of the mass in sinks versus time in 
fig- El [top]) is large from the start, and we observe no long 
period of slow, a ccelerating star format i on, co ntrary to 
the suggestion of iPalla fc Stahler I (|2000l 12002ft that star 
formation accelerates in time. Problems w i th th is sug- 
gestion have been discussed bv lHartmann I |2003), and, 
within the framework of their assumptions and limita- 
tions, our simulations do not confirm it. An important 
question is whether this result will persist when mag- 
netic fields are included. It is possible that during the 
early stages of the cloud's evolution, it may behave un- 
der magnetically subcritical conditions, giving very low 
star formation rates (SFRs), and then transit into a su- 
percritical regime, with higher SFRs f £)5.1.5[) . 

5.1.4. Latency period 

The clouds in our simulations do spend a long la- 
tency period between the beginning of the compres- 
sive motion that creates the cloud and the time at 
which they begin forming stars (~ 14-17 Myr in the 
three runs we have considered). However, this long 
period is most probably spent in the atomic phase, 
since the column densities of the star-forming re- 
gions in the simulations are co mparable to those re- 
quired for molecule formation ([Franco fc Cox I 119861 : 
lHartmann. Ballesteros-Paredes fc Bergin 1 120011 ). More 
specifically. iBergin et al. I (|2004ft have found that the 
timescale for CO molecule formation is essentially that 
required for reaching a dust extinction of Ay ~ 
0.7. Note also that during this time, the for- 
mation of H2 can be fast within density peaks 
(cores), even if most of the mass of the cloud is 
still in the low-density r egime (jPavlovski et al. 1 120021 : 
iGlover fc Mac Low Il2006l ). Thus, the concept of "rapid" 
star formation refers to the time elapsed between the ap- 
pearance of a (CO) molecular cloud and the onset star 



formation, as well as to the rate of the star formation 
process itself. 

Note that the latency period is in fact by definition of 
the order of the crossing time of the entire large-scale 
compr essive wave, in ag reement with observational evi- 
dence (|Elmegreen l l2000h . 

5.1.5. Implications for HI envelopes and magnetic criticality 

The formation of GMcs by an accumulation process 
such as the one modeled here has two important ad- 
ditional implications. First, because the process starts 
with lower-density gas (the WNM in our simulations) 
that is compressed by some external agent, GMCs, 
which are the "tip of the iceberg" of the density distri- 
bution, are expected to have in general HI envelopes, 
which would constitute the corresponding "body of the 
iceberg" . This gas would include transition material 
traversing the unstable phase between the warm to the 
cold phases, w hich would appear to have a nearly iso- 
baric behavior (|Gazql. Vazquez-Semadeni fc Kiml 20051: 



Audit fc Hennebelle I 2005; IVazquez-Semadeni et all 
2006ft. in agreement w i th th e conclusion reached by 
Andersson fc Wannier I (|1993ft from observations of 



HI envelopes of MCs. Such envelop es are routinely 
observed (e.g., iBlitz fc Thaddeus I fl980t IWanniei 
Lichten fc Morris]: lAndersson. Wannier fc Morris 

I1993F 



199lT 



Andersson & Wannier! 



Ballesteros-Paredes, Hartmann & Vazquez-Semadeni 
19991 iBrunt I l2003t IBlitz et al. I 120061 see also the 
combined HI and CO maps of Blitz & Rosolowsky 
2004). Weaker compressions than we have modeled 
may produce mostly thin CN M sheets, with little or no 
molecular g discussed in lVazquez-Semadeni et al~l 
(l20M ). 

Second, this scenario of molecular cloud formation im- 
plies that the mass-to-fiux ratio of the cloud is a vari- 
able quantity as the cloud evolves. This ratio is equiv- 
alent to the ratio of column density to magnetic field 
strength ( Na kano fc Nakamura I [19781 ) . with the critical 
colum density given by £ ~ 1.5 x 10 21 [B/5fiG] cm~ 2 . 
Although in principle under ideal MHD conditions the 
criticality of a magnetic flux tube involves all of the 
mass contained within it, in practice it is only the 
mass in the dense gas phase that matters, because 
the diffuse gas is not significantly self-gravitating at 
th e size scales of MC complexes. As pointed out 
bv lHartmann. Ballesteros-Paredes fc Bergin 1 (|2001l ). the 
above value of the dense gas' column density for magnetic 
criticality is very close to that required for molecule for- 
mation (eq. [T]) and for gravitational binding (eq. [2]), 
and therefore, the cloud is expected to become magneti- 
cally supercritical nearly at the same time it is becoming 
molecular and self-gravitating. 

5.2. Limitations and future work 
Our simulations are limited in a number of aspects. 



First and foremost, as already mentioned in §3.2| one 
major shortcoming of our simulations is the absence of 
feedback from the stellar objects onto the cloud. This is 
certainly an unrealistic feature. Second, our clouds are 
in a regime of turbulence decay after the inflows have 
subsided. This may or may not be an unrealistic feature, 
and rather it may represent a fraction of the population 
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of clouds. Finally, we have neglected magnetic fields al- 
together, due to the non-existence of suitable SPH algo- 
rithms including MHD for the problem of fully developed 
turbulent flows. 

All of these limitations tend to exaggerate the SFE in 
our simulations, as the turbulence within them is not re- 
plenished by either a continuing flow or by stellar en- 
ergy feedback. It is possible that a continued inflow 
(and therefore a sustained turbulence level) could pre- 
vent global collapse altogether, with only local collapse 
events happening randomly, although certainly the in- 
flows cannot last indefinitely. Thus, it would appear 
that global collapse is inescapable. However, once star 
formation begins, the stellar energy feedback is likely 
to either be able to halt the global collapse and dis- 
perse the cloud dFranco. Shore fe Tenorio-Tagle I Il994t 



lHartmann. Ballesteros-Paredes fc Bergin 2001) or else 
maintain it in rough equilibrium (e.g.. [ Matzner 1 1 2002 1: 
Nakamura fc Li 1 120051: iTan. Krumholz fe McKeel l200a 
Krumholz, Matzner fe McKee 1 12006ft . Thus, two im- 
portant questions to address in future papers are, one, 
what is the evolution like when star formation begins 
before the inflows subside. It is likely that this case 
will have much smaller SFEs. Two, whether the stel- 
lar energy feedback tends to disperse the clouds, or 
else to maintain them in rough equilibrium. Obser- 
vationally, star clusters older than 5-10 Myr tend to 
be al re ady devoid of gas dLeisawitz. Bash fc Thaddeus I 
19891 lHartmann. B allcstcros-Parcd es fc Bergin I 120011 : 



B allcsteros-Paredes fc Hartmann 



2006) , suggesting that 



the effect is more disruptive than equilibrating. 

Finally, our simulations have neglected the magnetic 
field. As discussed in 32J this is not a crucial omission 
if GMCs are in general magnetically supercritical, since 
the supercritical case is qualitatively equivalent to the 
nonmagnetic case, with the only difference being that it 
behaves as if it were less massive than a non- magneti c 
cloud of the same mass (|Shu I Il992t lHartmann! Il998f ). 
Thus, we consider that our models are still representa- 
tive of the large-scale evolution up to the early times 
of star formation in supercritical clouds. The main dif- 
ference is expected to be that magnetized, supercritical 
clouds should have lower SFEs than nonmagnet i c ones 
(IVazauez-Semadeni. Kim fc Ballesteros-Paredes I 120051 : 
iNakamura fc Lill2005l) . ' On the other hand, the evolu- 
tion of subcritical clouds will certainly differ from the 
models presented here, because they will have no global 
tendency to collapse. At any rate, it is necessary to per- 
form simulations of the full process in the presence of 
magnetic fields. This will presumably require the usage 
of AMR techniques incorporating the analogue of sink 
particles and stellar feedback, and will be pursued in fu- 
ture papers. 

6. SUMMARY AND CONCLUSIONS 

In this paper we have presented a suite of numeri- 
cal simulations designed to investigate jointly the for- 
mation of molecular clouds (MCs) and of stars within 
them. The simulations use an SPH scheme including 
self-gravity, sink particles and cooling leading to ther- 
mal bistability. Magnetic fields are neglected. The sim- 
ulations describe the collision of oppositely-directed gas 
streams ( "inflows" ) in the warm neutral medium (WNM) 
at moderately supersonic velocities (each with a velocity 



of 9.2 km s~ , or a Mach number of 1.22 in the unper- 
turbed WNM). Three simulations were considered, vary- 
ing the length of the inflow and the amplitude of the 
initial velocity fluctuations in the gas. The collisions trig- 
ger a transition to the cold phase in the gas, and simul- 
taneously generate turbulence in the resulting "cloud" , 
defined as the gas at densities n > 50 cm -3 . The in- 
flows secularly weaken in time, and so does the turbu- 
lence level in the cloud, implying that the turbulence 
in the clouds gradually transits from being continuously 
driven to being in a decaying regime. The cooling and 
the mass gain of the dense gas eventually cause the cloud 
to contain a large number of Jeans masses at the mean 
conditions. Moreover, because the cloud is supersonically 
turbulent, locally the density can be much larger, with 
a correspondingly lower Jeans mass and shorter free-fall 
times than those of the whole cloud. 

Thus, by the time the inflows have almost subsided, 
the cloud engages in global gravitational collapse, but 
shortly thereafter it begins to produce numerous local 
collapse events that occur on much smaller timescales 
because of the larger densities, so that by the time the 
global collapse of the cloud is completed, it has converted 
most of its mass to collapsed objects (sink particles). The 
sinks in general represent star clusters, due to the limited 
mass resolution an d lack of modeling of opacity-l imited 
fragmentation fe.g. lBate. Bonnell fc Bromm~l l2003) . Our 
simulations can only be considered reliable up to the time 
when the mass in collapsed objects (sinks) implies a high 
enough number of massive stars that they would disrupt 
the cloud. Nevertheless, the evolution up to that point 
shows a number of relevant results: 

1. The clouds are never in virial equilibrium during this 
period. Instead, they continually evolve, increasing their 
mass and self-gravitating energy, until the latter becomes 
comparable or larger than the thermal and kinetic ener- 
gies. Some ~ 5 Myr after global contraction began, star 
formation begins. 

2. In spite of not being in equilibrium, the near- 
equipartition condition \E S \ ~ 2Ey i is approximately sat- 
isfied as soon as \E g \ becomes comparable to Ek, because 
of the gravitational contraction of the cloud, with both 
quantities increasing simultaneously. This occurs long 
before the onset of star formation. This fact can explain 
the observed state of apparent virialization of GMCs. 

3. The near equipartition is a signature of global gravi- 
tational collapse, not equilibrium, and suggests a return 
to Goldreich & Kwan's (1974) scenario of global grav- 
ita tional contraction in MCs. However, the criticism 
bv lZuckerman fc Palmer 1 ([T974), namely that an exces- 
sively large star formation should result through this pro- 
cess, is avoided in part because the nonlinear turbulent 
density fluctuations collapse earlier than the whole cloud, 
involving only a fraction of the total mass, and in part be- 
cause as soon as the stars form they probably contribute 
to dispersing the cloud, or at least halting its global col- 
lapse. Further reduction of the SFE may occur in the 
presence of supercritical magnetic fields. 

4. Local collapse events begin to occur after global col- 
lapse has begun, but they occur rapidly, requiring only 
~ 3 Myr to convert ~ 15% of the cloud's mass into stars. 
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5. The star formation rate, measured by the slope of the 
mass in stars versus time, is large from the beginning. 
Within the framework, assumptions and limitations of 
our simulations, no long period of slow, accelerating star 
formation is observed. 

6. Longer inflow durations maintain larger turbulent ve- 
locity dispersions in the clouds and delay the onset of 
both global and local collapse. Instead, larger amplitudes 
of the initial velocity fluctuations have little effect in de- 
laying the collapse. The latter effect can be attributed 
to the fact that the turbulence is already dissipating in 
the clouds by the time global collapse begins. 

7. A long period (~ 14-17 Myr) of "dormancy" does oc- 
cur between the time when the cloud begins to form and 
the time when star formation begins. Nevertheless, it is 
likely that most of this time is spent with the gas being 
in atomic form, since the column density of star form- 
ing regions in our simulations are comparable to values 
required for molecular gas formation, as reported in the 
literature. Thus, our simulations support the notion that 
star formation occurs almost simultaneously with the for- 
mation of molecular gas, under local ISM conditions. 

We conclude that our simulations support the 
scenario of rapid star formation after molecu- 
lar gas has formed, involving accumulations of 
gas from distances of a few hundred parsecs 
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TABLE 1 

Run parameters. 



Run name 


ibox a (pc) 


4nf b (pc) 


v ini c (km s 1 ) 


-Mi„f d 


frms,i C (km S 1 ) 


M bo J (M Q ) 




AM 11 (M©) 


L128A«0.24 


128 


48 


9.20 


1.22 


0.24 


6.582 x 10 4 


1.643 x 10 6 


0.04 


L128At;0.66 


128 


48 


9.20 


1.22 


0.66 


6.582 x 10 4 


1.643 x 10 6 


0.04 


L256Au0.17 


256 


112 


9.20 


1.22 


0.17 


5.253 x 10 5 


3.242 x 10 6 


0.16 



a Physical size of computational domain.* 3 Linear size of each inflow. c Speed of inflows. d Mach number of inflows, with respect to sound 
speed of unperturbed WNM. e One-dimensional rms speed of initial imposed velocity fluctuations/ Total mass in the computational 
domain. 8 Number of SPH particles. 11 Mass resolution, equal to the mass of each SPH particle. 



TABLE 2 

Column densities of star forming regions in run L256A^0.17. 



Region Time (Myr) Position in yz plane (pc) N ( cm 2 ) 



1 17.0 (144,135) 

2 17.26 (131,142) 

3 18.06 (112,125) 

4 18.59 (137,115) 



9.58 x 10 20 
2.27 x 10 21 
4.89 x 10 20 
1.52 x 10 21 
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Fig. 1. — Setup of the simulations. Within a cubic box measuring Lbox on a side, two cylindrical inflows, each of length £] n f, radius 
r; n f, and aligned with the x-axis of the domain, are set to collide. The inflows have a velocity corresponding to -Mi n f = 1.22 with respect 
to the undisturbed medium. Also, a fluctuating velocity field with amplitude t> r ms,i is added at the begining of the run in order to trigger 
the dynamical instability of the resulting compressed layer. Table [T] lists the values used for these parameters. 
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Fig. 2. — Thermal-equilibrium pressure versus density for the cooling and heating functions given by eqs. <[3j and J4}. The vertical dotted 
line indicates the uniform initial density in the simulations. The tangent to the curve at that location has a near-zero slope, indicating 
7c Ri 0. At n = 0.2 cm -3 , 7 C = 0.83. 




E 



50 100 150 200 250 50 100 150 200 



Fig. 3. — Edge-on views of run L256Au0.17 at t = 2.66 Myr (left) and 20.0 Myr (right). The arrows indicate the velocity field. This 
figure is available as an mpeg animation in the electronic edition of The Astrophysical Journal. 




Fig. 4. — Face-on views of run L256A?j0.17 at t = 17.0 Myr (left) and 18.9 Myr (right). The dots indicate the projected positions of the 
sink particles This figure is available as an mpeg animation in the electronic edition of The Astrophysical Journal.. 
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Fig. 5. — Top: Evolution of the cloud, sink, and total (cloud+sink) mass for run L256At)0.17. Bottom: Evolution of the velocity 
dispersion of the dense gas (n > 50 cm -3 ) in each of the three coordinate axes. From t ~ 6 Myr through fa 12 Myr, a x decreases, reflecting 
the smooth end of the inflows. At this time, the gravitational collapse of the cloud (in the yz plane) sets in, generating an increase in cr y 
and o~ z ■ While initially unaffected by this collapse, a x increases starting at t fs 18 Myr, coinciding with the onset of star formation. 
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Fig. 6. — Evolution of th e cloud mass and of the minimum mass that is not dispersed by stellar ionizing radiation under the two estimates 
given by eqs. J9} and dlOI I for run L256A?;0.17. The cloud is expected to be disrupted when the minimum cloud mass for non-disruption 
exceeds the actual cloud mass, at the time indicated by the vertical dotted line. 
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Fig. 7. — Evolution of the mass- weighted density {left frame) and temperature (right frame) histograms (normalized to the maximum) 
for run L256A«0.17 at various times. 
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Fig. 8. — Top: Evolution of the total gravitational energy (in absolute value) of the numerical box and of the kinetic and thermal energies 
in a cylinder of radius 32 pc and length 16 pc centered in the middle of the numerical box for run L256Ai;0.17. Bottom: Similar to the 
top panel, but with the thermal and kinetic energies calculated for the dense gas (n > 50 cm -3 ) only. 
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Fig. 9. — Evolution of the velocity dispersion in the dense gas (n > 50 cm 3 ) along each of the three coordinate axes for run L128A^0.24. 
Compare to fig. \5\ (right), noting the different extension of the axes. 



